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ABSTRACT 

Recently,  Energetic  Ionic  Liquids  (EILs)  have  emerged  as  potential  alternative  hypergolic 
propellants  to  replace  toxic  monomethyl  hydrazine  (MMH).  Ionic  liquids  have  no  appreciable 
vapor  pressures,  are  safe  to  handle  and  can  be  tailored  to  have  desired  properties.  One 
challenge  is  to  design  EILs  with  short  ignition  delays.  In  this  paper,  we  report  density  functional 
tight  binding  molecular  dynamics  simulations  (DFTB-MD)  to  study  initial  stages  of  hypergolic 
reaction  between  an  EIL  and  nitric  acid.  Calculations  were  performed  at  various  temperatures, 
and  reaction  mechanisms  were  identified.  The  reaction  products,  such  as  H20,  HNCO  and  C02, 
predicted  by  DFTB-MD  simulations  are  in  agreement  with  the  recent  experiments  by  AFRL. 
These  simulations  predict  that  the  reaction  mechanism  is  very  complex  and  it  changes  with 
temperature. 


INTRODUCTION 

Hypergolic  propellants  have  been  used  in  small  to  medium  range  rocket  engines  since  World  War 
II.  This  category  of  liquid  propellant  ignites  almost  instantaneously  when  exposed  to  powerful 
oxidizers,  such  as  IFRNA.  Currently  the  most  widely  used  hypergolic  propellant  is  MMH. 
Although  monomethyl  hydrazine  (MMH)  has  been  used  for  decades,  it  is  highly  toxic  and 
carcinogenic.  It  also  has  high  vapor  pressure  at  room  temperature  which  makes  it  difficult  to 
handle  and  transport.  In  addition,  recent  environment  and  health  concerns  encourage  the 
development  of  environment  friendly  fuels.  Therefore,  one  of  the  major  goals  of  the  Air  Force,  and 
in  general,  DoD  is  to  search  for  hypergolic  propellants  that  can  replace  MMH  without  sacrificing 
its  performance.  Recently  Energetic  Ionic  Liquids  (EIL)  have  emerged  as  an  attractive  alternative 
due  to  their  potential  applications  as  hypergolic  bipropellants1'5.  Due  to  their  negligible  vapor 
pressure,  they  are  more  environment-friendly  and  safe  to  handle  and  transport  than  the  currently 
used  hydrazine  based  fuels.  In  addition,  EILs  are  amenable  to  the  modification  of  their  physical 
properties  such  as  density  and  melting  point  via  the  introduction  of  different  chemical 
functionalities.  For  these  reasons,  EILs  are  being  seriously  considered  as  next  generation 
hypergolic  propellants. 

The  most  important  aspect  of  a  hypergolic  fuel  is  its  ignition  delay  (ID)  time  when  mixed  with  an 
oxidizer,  such  as  IRFNA.  The  longer  the  ID  is,  the  larger  the  combustion  chamber  must  be  to 
avoid  pressure  spikes  that  could  rupture  the  engine.  So,  a  short  ID  is  a  necessary  condition  for 
hypergolic  fuel.  Therefore,  the  current  challenge  is  to  design  a  hypergolic  fuel  which  is  high- 
performing,  nontoxic,  safe  to  handle  and  transport,  has  ID  comparable  to  MMH  (~5  ms)  and  can 
be  used  with  existing  propulsion  systems  with  no  or  minor  modifications.  Although  EILs  have 
come  long  way  towards  meeting  these  requirements,  their  ID  is  still  considerably  longer  than 
MMH.  Therefore,  one  challenge  is  to  design  EILs  with  short  ignition  delay.  Currently  there  is  no 
technique  to  predict  ID  prior  to  synthesis  of  EILs.  Very  recently,  we  have  developed  a  preliminary 
model  for  ID  prediction  using  a  Quantitative  Structure  Property  Relationship  (QSPR)  and  Artificial 
Neural  Network  (ANN).  Although  these  predictive  models  are  useful  for  design  purposes,  it  does 
not  provide  fundamental  understanding  of  the  reaction  mechanism. 

In  order  to  obtain  a  fundamental  understanding  of  these  reaction  steps,  quantum  based  methods 
are  the  methods  of  choice.  The  quantum  based  methods  are  quite  complex,  time  consuming  and 
often  require  multi-disciplinary  expertise.  For  example,  ab  initio  quantum  chemical  calculations 
when  combined  with  reaction  rate  theory  and  kinetics  modeling  can,  in  principle,  be  capable  of 
predicting  the  ID.  However,  such  a  procedure  is  extremely  time  consuming  (may  take  years  for 
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each  EIL),  and  each  EIL  must  be  studied  on  a  case  by  case  basis.  Therefore,  such  a  method  is 
not  practical,  but  can  be  useful  in  understanding  some  of  the  fundamental  aspects  of  ignition 
reactions.  Quantum  molecular  dynamics  (MD),  when  used  in  conjunction  with  semi-empirical 
forcefields,  is  faster  than  ab  initio  calculations,  and  reveals  the  dynamics  of  bond-making  and 
breaking.  In  addition,  such  calculations  can  be  done  for  several  hundred  of  atoms  with  minimal 
computer  resources,  and  a  realistic  nitric  acid  environment  can  be  set  around  an  EIL  molecule  for 
accurate  predictions.  This  is  in  contrast  to  high  level  ab  initio  methods  where  calculations  are 
inherently  with  vapor  phase  approximations.  Although,  at  present,  quantum  reactive  MD 
calculations  cannot  be  performed  at  the  millisecond  time  scale  (i.e.  the  time  scale  of  ID),  one  can, 
at  least,  obtain  some  insight  into  the  initial  stages  of  hypergolicity.  We  therefore  preferred 
quantum  reactive  MD  simulation  with  density  functional  tight-binding  (DFTB)  forces  over 
traditional  ab  initio  methods  to  obtain  fundamental  understanding  of  hypergolic  reactions.  We  will 
refer  this  technique  as  DFTB-MD  in  our  following  discussion. 

In  this  paper,  we  first  outline  the  DFTB-MD  method.  We  then  perform  DFTB-MD  simulations  at 
various  temperatures,  identify  the  predicted  reaction  products  and  reaction  pathways  at  various 
temperatures. 


RESULTS  AND  DISCUSSION 


Computational  Method 

One  of  the  major  advantages  of  the  this  method  is  that  the  forces  required  for  atom  movements 
are  computed  using  a  self-consistent  charge  tight  binding  method  (in  the  framework  of  density 
functional  theory)  on  the  fly  without  using  any  empirical  analytical  expressions.  Since  it  is  based 
on  quantum  formalism,  the  chemical  reactions  are  naturally  taken  into  account,  and  no  special 
treatment  is  required  for  chemical  reactions.  Recently,  the  same  technique  has  been  successfully 
applied  to  study  the  dynamics  of  HMX  decomposition  in  extreme  conditions6.  To  our  knowledge, 
this  is  the  first  report  of  applying  DFTB-MD  to  study  hypergolicity  between  an  EIL  and  nitric  acid. 
In  the  following,  we  outline  the  main  features  of  the  DFTB  method,  and  we  refer  to  the  literature 
for  additional  details7'10. 

■s  The  DFTB  method  is  based  on  the  second  order  expansion  of  the  Kohn-Sham  energy 
functional  from  DFT  with  respect  to  charge  density  fluctuation  relative  to  a  reference 
density. 

S  Unlike  the  traditional  tight-binding  method,  the  DFTB  calculates  the  charge  transfer 
between  atoms  and  performs  population  analysis  similar  to  that  done  in  other 
sophisticated  quantum  chemical  methods. 

■S  Like  DFT,  DFTB  uses  non-orthogonal  atomic  orbitals. 

•s  The  repulsive  part  of  the  energy  functional,  which  is  a  sum  of  inter-atomic  two-body 
potentials,  is  obtained  by  fitting  with  small  molecules  using  DFT  data. 

MD  simulations  were  performed  with  constant  volume  and  a  fixed  number  of  atoms.  Two  types  of 
simulations  were  performed:  1)  increase  in  temperature  from  10K  to  3000K  with  time,  and  2)  at  a 
constant  temperature.  The  increasing  temperature  simulation  was  performed  in  order  to  examine 
the  sequence  of  events  that  happens  as  the  temperature  changes  rapidly.  The  time  step  used  for 
MD  simulation  was  chosen  to  be  1  fs.  At  each  time  step,  DFTB  forces  were  calculated,  and 
Velocity  Verlet  algorithm  was  used  to  determine  the  position  of  the  atoms  in  the  following  time 
step.  All  simulations  were  run  at  least  up  to  10  ps. 

Model  Construction 


We  chose  1-allyl-3methyl  imidazolium  dicyanamide,  an  EIL,  and  pure  HN03i  an  oxidizer,  to  study 
the  hypergolic  reactions.  In  order  to  study  the  reactions  correctly,  it  was  critical  to  construct  an 
appropriate  and  a  reasonably  realistic  model.  o  To  accomplish  this,  we  first  constructed  a 
simulation  box  with  dimension  of  24A  x  12A  x  12A  filled  with  nitric  acid.  An  appropriate  number  of 
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nitric  acid  molecules  were  put  into  the  box  so  that  it  approximately  reproduced  the  experimental 
density  of  the  nitric  acid  (1.51  gm/cc).  The  total  number  of  nitric  acid  molecules  therefore 
considered  was  50.  Periodic  boundary  condition  was  applied  so  that  when  an  atom  left  the  box 
during  the  simulation,  another  one  was  inserted  in  an  equivalent  position  from  the  opposite 
direction.  The  preliminary  nitric  acid  configuration  was  optimized  for  100  optimization  steps  in 
DFTB  in  order  to  obtain  a  relaxed  configuration.  However  this  resulted  in  an  ordered  structure, 
instead  of  a  disordered  liquid  structure.  To  obtain  a  liquid  structure,  we  gradually  heated  up  the 
system  from  OK  to  600K,  where  the  system  melted,  and  then  cooled  down  to  room  temperature 
(298.1 5K).  This  procedure  yielded  a  disordered  liquid  structure  for  nitric  acid,  as  shown  in  Figure 
1. 

We  then  manually  inserted  one  molecule  of  EIL  in  the  middle  of  the  simulation  box  by  removing 
three  molecules  of  nitric  acid.  Such  an  adjustment  kept  the  density  of  the  system  nearly 
unchanged.  The  entire  configuration,  except  the  atoms  in  EIL,  was  further  relaxed  in  order  to 
minimize  close  contacts.  The  atoms  in  EIL  were  fixed  to  prevent  any  bond  breaking  or  forming 
during  the  minimization  steps.  Before  the  EIL  was  inserted,  the  gas  phase  geometry  of  the  EIL 
was  optimized  with  DFTB  and  charges  were  examined.  The  charges  on  the  cation  and  anion 
were  found  to  be  very  close  to  +1  and  -1,  respectively  confirming  that  the  DFTB  charge 
distributions  were  correct.  The  EIL  containing  box  filled  with  nitric  acid  was  taken  as  starting  point 
of  our  MD  simulations  (Figure  2). 

It  should  be  noted  that  IDs  for  hypergolic  reactions  with  nitric  acid  are  of  the  order  of  10 "3s  and 
longer.  Reactive  MD  simulation  with  such  time  scale  has  not  been  reported  to  date,  and  it  is 
impractical  to  cover  such  a  time  scale  with  standard  MD  techniques.  However,  coupling  of  DFTB 
with  Temperature  Accelerated  Dynamics  (TAD)11  could  allow  simulation  of  such  a  time  scale.  It 
should  be  pointed  out  that  the  TAD  is  not  the  same  as  artificially  increasing  the  temperature  to 
enhance  the  reaction  rate,  which  is  used  for  the  reported  simulations  in  this  study  and  is  the  most 
widely  used  technique  to  study  slow  chemical  processes.  The  goal  of  the  present  investigation, 
therefore,  is  not  to  study  the  progress  of  hypergolic  reactions  over  the  real  time  scale,  but  to  get 
an  idea  of  the  reactions  and  intermediates  that  might  be  involved  during  the  initial  stages  of 
reaction. 

Simulation  Results 


As  mentioned  earlier,  two  different  types  of  simulations  were  performed:  1)  increasing  the 
temperature  with  time  and  2)  at  constant  temperature. 

Simulation  with  Increase  in  Temperature  with  Time  to  3000K: 

In  this  simulation,  the  temperature  of  the  system  was  increased  from  10K  to  3000K  over  a  time 
period  of  7  ps,  and  then  kept  constant  for  another  3  ps.  The  objective  was  to  follow  the  sequence 
event  as  the  temperature  increased. 

Figure  3  shows  the  reactions  that  occurred  during  the  time  of  the  simulation.  Figure  4  shows  the 
snapshots  at  different  times  when  these  reactions  occurred.  The  dicyanamide  (DCA)  anion 
picked  up  two  protons  and  was  simultaneously  oxidized,  then  split  into  HNCO  and  HNCN.  The 
formation  of  HNCO  occurred  at  around  7  ps  with  the  current  simulation  conditions.  On  the  other 
hand,  the  methyl  group  in  the  side  chain  was  oxidized  to  -CH2OH.  On  the  other  side  chain,  the 
terminal  double  bond  rearranged  (from  3  position  to  2  position).  The  HNCN,  which  formed  from 
the  DCA,  attacked  at  the  side  chain  of  the  cation.  This  was  followed  by  a  series  of  reactions 
(oxidation  and  bond  breaking)  leading  to  acetylene  and  CO  (via  HCO).  HCO  (which  is  the 
precursor  for  CO  formation)  was  formed  along  with  HNCN  which  appeared  to  be  responsible  for 
rupturing  the  cation  completely.  HNCN  finally  was  converted  to  H2N-CN  and  remained  stable 
until  end  of  simulation  (10  ps).  One  important  observation  in  this  simulation  was  the  formation  of 
HNCO  from  DCA  which  was  experimentally  observed  by  the  AFRL  group.2  We  have  not 
observed  any  C02  and  N20  formation  during  the  length  of  the  simulation.  However,  there  is  a 
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possibility  that  CO  could  further  be  oxidized  to  C02  at  a  much  longer  time  scale.  In  this 
simulation  we  did  not  observe  the  formation  of  oxides  of  nitrogen  from  the  IL,  but  their  formation 
from  nitric  acid  could  not  be  ruled  out  since  there  were  many  oxidation  reactions  in  the  system 
where  one  of  the  oxygens  of  HN03  was  donated.  In  all  MD  simulations  presented  here,  we  have 
identified  and  tracked  the  intermediates  and  products  by  visual  inspection  of  the  trajectories. 

Simulation  at  1000K: 

Simulation  was  performed  at  a  constant  temperature  of  1000K  with  the  same  initial  condition  as 
before.  The  DCA  anion  quickly  picked  up  a  proton  to  form  NC-NH-CN.  However,  anion  and 
protonated  DCA  did  not  react  during  the  time-scale  of  the  simulation  (1  Ops).  This  clearly 
indicates  that  the  time-scale  of  reaction  is  too  large  compared  to  that  of  MD,  and  reactions  could 
not  be  observed  unless  the  temperature  is  artificially  increased  or  the  TAD  methodology  is 
implemented  and  used. 

Simulation  at  2000K: 

At  2000K,  as  observed  for  1000K,  protonation  of  DCA  occurred.  The  unsaturated  side  chain  of 
the  cation  was  oxidized  and  detached  from  the  imidazolium  group.  The  products  formed  during 
the  time  of  simulations  were  an  unsaturated  aldehyde  and  1 -methyl  imidazolium  ion  (Figure  5). 

Simulation  at  2500K: 

Unlike  at  1000K  and  2000K,  many  events  were  observed  at  2500K  (Figure  6).  The  anion 
captured  a  proton  to  form  NC-NH-CN  and  then  HCH-N-CN.  The  protonated  anion  eventually  was 
oxidized  to  form  HNCO  and  HNCN  via  a  series  of  rearrangement  and  decomposition  reactions. 
The  final  step  for  HNCO  formation  happened  at  9.98  ps.  On  the  other  hand,  the  methyl  side  chain 
of  the  cation  was  oxidized  to  form  formaldehyde  (H2C=0)  which  led  to  CO  via  HCO  formation  (at 
2.92  ps).  The  five-membered  imidazole  compound  (remained  after  the  oxidation  of  the  methyl 
group)  underwent  a  series  of  reactions,  such  as  oxidation,  hydrogen  abstraction  and  ring  rupture. 
The  rupture  of  the  C-N  bond  in  the  ring  occurred  via  reaction  with  HNCN  which  formed  from  the 
anion.  We  continued  this  simulation  for  25  ps  in  order  to  observe  the  fate  of  the  five-membered 
ring.  The  most  important  finding  was  the  formation  of  C02  via  oxidation  of  the  imidazole 
compound.  We  also  observed  some  H20  as  byproduct.  The  important  chemical  transformations 
are  shown  in  Figure  7. 

Simulation  at  3000K: 

Reactions  observed  at  3000K  simulation  (constant  temperature)  are  shown  in  Figure  8.  Figure  9 
shows  the  important  structures  formed  at  different  times  during  the  simulation.  At  3000K,  the  DCA 
anion  captured  a  proton  (NC-NH-CN  ->  NC-N-CNH),  and  then  abstracted  a  hydrogen  from  H20 
(formed  via  2HN03  -->  N02  +  N03  +  H20)  to  form  a  short-lived  HNC-N-CNH  species  which 
immediately  formed  HNC  and  HNCN.  HNC  finally  was  oxidized  to  NO  and  CO  via  a  series  of 
reactions  including  HNCO.  In  fact,  HNCO  formed  at  3.56  ps,  lived  until  5.52  ps  and  then  was 
converted  to  NCO.  NO  and  CO  formation  occurred  at  6.74  ps  at  this  temperature.  One  other 
hand,  allyl  side  chain  of  the  cation  lost  two  successive  hydrogens  to  form  a  triple  bond.  The 
methyl  side  chain  was  oxidized  to  formaldehyde  which  eventually  produced  CO  (H2C=0  ->  HCO 
->  CO).  Formation  of  CO  from  cation  occurred  much  faster  (0.92  ps)  than  that  from  the  anion. 
The  cyclic  part,  present  after  oxidation  of  the  methyl  side  chain,  underwent  several  reactions:  a) 
ring  opening,  b)  ring  expansion  to  form  a  pyridine-like  structure,  c)  rupture  of  the  pyridine-like 
structure  via  oxidation  and  d)  finally  the  formation  of  CO,  NO  and  HCN.  A  CN  fragment, 
generated  from  the  pyridine-like  structure,  reacted  with  HNCN  (formed  from  the  anion)  to  form 
HNC-C-CN  which  eventually  was  converted  to  HNCO  via  oxidation  at  9.72  ps.  The  major  finding 
of  this  simulation  is  the  mechanism  for  formation  of  HNCO,  NO  and  CO.  HNCO  and  CO  formed 
from  both  anion  and  cation  of  the  IL.  Also,  like  at  2500K,  we  observed  H20  formation.  One  of  the 
most  interesting  aspects  of  these  simulations  is  that  the  mechanisms  of  HNCO,  CO  and  NO 
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formation  were  different  at  different  temperatures.  This  gives  some  indication  that  the  hypergolic 
reaction  of  the  EIL  could  be  more  complex  than  thought  before. 

SUMMARY  AND  CONCLUSIONS 

This  paper  reports,  for  the  first  time,  DFTB-MD  simulations  to  understand  the  initial  stages  of 
hypergolic  reactions  between  an  EIL  and  nitric  acid.  In  these  simulations,  an  EIL  molecule  was 
placed  in  a  bath  of  nitric  acid.  The  simulations  were  performed  at  various  temperatures  (such  as 
1000K,  2000K,  2500K  and  3000K),  and  the  products  and  reactions  were  identified.  The  products 
identified  were  HNCO,  CO,  NO,  N02,  C02  and  H20.  It  should  be  noted  that  HNCO  and  C02  were 
observed  by  the  AFRL  group  in  their  experimental  study.2  However,  we  did  not  observe  N20  (as 
observed  by  experiments)  during  the  timescale  of  the  simulations.  We  anticipate  that  N20  will 
form  at  much  later  stages  in  the  process.  DFTB-MD  simulations  show  the  reaction  between  the 
EIL  and  nitric  acid  to  be  very  complex,  and  that  the  mechanism  appears  to  be  significantly 
different  at  different  temperatures. 
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FIGURES 


Figure  1 :  Liquid  Structure  of  Nitric  Acid  Figure  2: 1-allyl  3methyl  dicyanamide  EIL  i 

HN03.  This  is  the  Starting  Configuration  for 
Simulations. 
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Figure  3:  Reactions  Observed  During  Simulations  with  Gradual  Increase  in  Temperature  to 
3000K.  Intermediates  Shown  Above  are  Actually  Observed  During  the  Simulation. 
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Figure  5:  Oxidation  Reaction  Occurred  at  2000K.  Intermediates  Shown  Above  are  Actually 

Observed  During  the  Simulation. 
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Figure  6:  Reactions  Observed  During  MD  Simulation  at  2500K.  Intermediates  Shown  Above  are 

Actually  Observed  During  the  Simulation. 
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Figure  7:  Snapshots  at  Different  Times  for  Important  Structures  Shown  in  Figure  6.  They  are 
Arranged  According  to  the  Progress  in  Time. 
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Figure  8:  Reactions  Observed  during  MD  Simulation  at  3000K.  Intermediates  Shown  Above  are  Actually 

Observed  During  Simulation. 
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Figure  9:  Snapshots  at  Different  Times  for  Selected  Important  Structures  Shown  in  Figure  8.  They 
are  Arranged  According  to  the  Progress  in  Time. 
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